rm(list = ls())
options(warn=-1)
library(rgl)
library(misc3d)
knitr::knit_hooks$set(webgl = hook_webgl)

path.bin <- './molecular-atlas-master/bin'
path.matrices <-'./molecular-atlas-master/data'

load(paste(path.matrices ,'atlasspots.RData',sep='/'))
load(paste(path.matrices , 'vivid-colors.RData', sep='/'))

if(!exists('aligned.atlas'))
  load(paste(path.matrices ,'alignedAtlas.RData',sep='/'))

if(!exists('atlas.stereo'))
  load(paste(path.matrices ,'atlasStereo.RData',sep='/'))

source(paste(path.bin,'plotFunctions.R',sep='/'))
source(paste(path.bin,'execFunctions.R',sep='/'))
source(paste(path.bin,'araAtlasFunctions.R',sep='/'))
source(paste(path.bin,'smoothingFunctions.R',sep='/'))
source(paste(path.bin,'layerFunctions.R',sep='/'))
source(paste(path.bin,'icFunctions.R',sep='/'))
source(paste(path.bin,'tsneFunctions.R',sep='/'))
source(paste(path.bin,'similarityIndexFunctions.R',sep='/'))
source(paste(path.bin,'scRemapFunctions.R',sep='/'))
source(paste(path.bin,'constantParameters.R',sep='/'))
source(paste(path.bin,'allenAnnotationsFunctions.R',sep='/'))
source(paste(path.bin,'meshCutsFunctions.R',sep='/'))

Load results and 3D locations of spots.

#STitch3D's results and CCFv3 coordinates.
load("spotstable.RData")
load("VOLUMESMALL.RData")
load("VOLUME.RData")
for (i in (0:34)){
  prop <- read.table(paste0("./res_mouse_brain/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
library(stringr)
prop_all$X <- str_split_fixed(prop_all$X, "-", 2)[, 1]
spots.table <- spots.table[prop_all$X, ]
spots.table$X <- rownames(spots.table)
spots.table <- merge(spots.table, prop_all, by=c("X"))
spots.table.model <- spots.table

#STitch3D's results and constructed 3D locations of spots.
load("spotstable.RData")
for (i in (0:34)){
  prop <- read.table(paste0("./res_mouse_brain/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
prop_all$X <- str_split_fixed(prop_all$X, "-", 2)[, 1]
spots.table <- spots.table[prop_all$X, ]
spots.table$X <- rownames(spots.table)
spots.table <- merge(spots.table, prop_all, by=c("X"))
colnames(spots.table)[10] <- "x_pixel"
colnames(spots.table)[11] <- "y_pixel"
coor_3d <- read.table(paste0("./res_mouse_brain/3D_coordinates.csv"), sep=",", header=TRUE)
coor_3d$X <- str_split_fixed(coor_3d$X, "-", 2)[, 1]
spots.table <- merge(spots.table, coor_3d, by=c("X"))
spots.table.ours <- spots.table

hpc_colors = c("#d62728", "#FFFF00", "#1f77b4", "#2ca02c")

STitch3D’s cell-type deconvolution result reconstructs curve-shaped distributions of four hippocampal neuron types in three cornu ammonis areas and dentate gyrus.

spots.table <- spots.table.ours

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   1

par3d(persp)

## NULL

um <- c(0.98223096, -0.05770408, -0.1785820, 0,
        0.02726252, -0.89759272, 0.4399813, 0,
        -0.18568322, -0.43703181, -0.8800704, 0,
        0.0000000, 0.0000000, 0.0000000, 1) # view1
view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))

idx_CA1 = (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA1>0.2)
idx_CA2 = (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA2>0.2)
idx_CA3 = (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA3>0.2)
idx_DG1 = (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_DG1>0.2)
spheres3d(spots.table[idx_CA1,]$x, spots.table[idx_CA1,]$y, spots.table[idx_CA1,]$z,
          col = hpc_colors[1], radius=0.5, 
          alpha=spots.table[idx_CA1,]$Ext_Hpc_CA1)
spheres3d(spots.table[idx_CA2,]$x, spots.table[idx_CA2,]$y, spots.table[idx_CA2,]$z,
          col = hpc_colors[2], radius=0.5, 
          alpha=spots.table[idx_CA2,]$Ext_Hpc_CA2)
spheres3d(spots.table[idx_CA3,]$x, spots.table[idx_CA3,]$y, spots.table[idx_CA3,]$z,
          col = hpc_colors[3], radius=0.5, 
          alpha=spots.table[idx_CA3,]$Ext_Hpc_CA3)
spheres3d(spots.table[idx_DG1,]$x, spots.table[idx_DG1,]$y, spots.table[idx_DG1,]$z,
          col = hpc_colors[4], radius=0.5, 
          alpha=spots.table[idx_DG1,]$Ext_Hpc_DG1)
spheres3d(spots.table$x, spots.table$y, spots.table$z,
          col = 'gray', radius=0.5, 
          alpha=0.02)
decorate3d()

spots.table <- spots.table.model

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   3

par3d(persp)

## NULL

um <- c(-0.6389505, -0.1887468, -0.7457324, 0,
        0.4022625, -0.9083022, -0.1147684, 0,
        -0.6556883, -0.3733116, 0.6562859, 0,
        0, 0, 0, 1) # optimal view1
view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
drawScene.rgl(list(VOLUMESMALL))
smp.AP <- 0
smp.ML <- 0
idx_CA1 = (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_CA1>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA1>0.2)
idx_CA2 = (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_CA2>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA2>0.2)
idx_CA3 = (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_CA3>spots.table$Ext_Hpc_DG1) &
    (spots.table$Ext_Hpc_CA3>0.2)
idx_DG1 = (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA1) &
    (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA2) &
    (spots.table$Ext_Hpc_DG1>spots.table$Ext_Hpc_CA3) &
    (spots.table$Ext_Hpc_DG1>0.2)
# spheres3d(spots.table[idx_CA1, ]$AP.paxTOallen - 530/2 + smp.AP, 
#           -spots.table[idx_CA1, ]$DV * 1000/25 - 320/2, 
#           spots.table[idx_CA1, ]$ML * 1000/25 + smp.ML, 
#           col = hpc_colors[1], radius=5, 
#           alpha=spots.table[idx_CA1,]$Ext_Hpc_CA1)
spheres3d(spots.table[idx_CA1, ]$AP.paxTOallen - 530/2 + smp.AP, 
          -spots.table[idx_CA1, ]$DV * 1000/25 - 320/2, 
          spots.table[idx_CA1, ]$ML * 1000/25 + smp.ML, 
          col = hpc_colors[1], radius=5, 
          alpha=1)
# spheres3d(spots.table[idx_CA2, ]$AP.paxTOallen - 530/2 + smp.AP, 
#           -spots.table[idx_CA2, ]$DV * 1000/25 - 320/2, 
#           spots.table[idx_CA2, ]$ML * 1000/25 + smp.ML, 
#           col = hpc_colors[2], radius=5, 
#           alpha=spots.table[idx_CA2,]$Ext_Hpc_CA2)
spheres3d(spots.table[idx_CA2, ]$AP.paxTOallen - 530/2 + smp.AP, 
          -spots.table[idx_CA2, ]$DV * 1000/25 - 320/2, 
          spots.table[idx_CA2, ]$ML * 1000/25 + smp.ML, 
          col = hpc_colors[2], radius=5, 
          alpha=1)
# spheres3d(spots.table[idx_CA3, ]$AP.paxTOallen - 530/2 + smp.AP, 
#           -spots.table[idx_CA3, ]$DV * 1000/25 - 320/2, 
#           spots.table[idx_CA3, ]$ML * 1000/25 + smp.ML, 
#           col = hpc_colors[3], radius=5, 
#           alpha=spots.table[idx_CA3,]$Ext_Hpc_CA3)
spheres3d(spots.table[idx_CA3, ]$AP.paxTOallen - 530/2 + smp.AP, 
          -spots.table[idx_CA3, ]$DV * 1000/25 - 320/2, 
          spots.table[idx_CA3, ]$ML * 1000/25 + smp.ML, 
          col = hpc_colors[3], radius=5, 
          alpha=1)
# spheres3d(spots.table[idx_DG1, ]$AP.paxTOallen - 530/2 + smp.AP, 
#           -spots.table[idx_DG1, ]$DV * 1000/25 - 320/2, 
#           spots.table[idx_DG1, ]$ML * 1000/25 + smp.ML, 
#           col = hpc_colors[4], radius=5, 
#           alpha=spots.table[idx_DG1,]$Ext_Hpc_DG1)
spheres3d(spots.table[idx_DG1, ]$AP.paxTOallen - 530/2 + smp.AP, 
          -spots.table[idx_DG1, ]$DV * 1000/25 - 320/2, 
          spots.table[idx_DG1, ]$ML * 1000/25 + smp.ML, 
          col = hpc_colors[4], radius=5, 
          alpha=1)